library(tidyverse)
library(ggpubr)
library(patchwork)

parse_qpcr_cp <- function(file, gene = NULL, group = NULL, transposed = FALSE){
  res_parse <- read_delim(file, skip = 1) |>
    select(Pos, Cp) |>
    filter(!is.na(Cp)) |>
    mutate(col_index = str_extract(Pos, '\\d+') |> as.numeric(),
           row_index = str_extract(Pos, '[A-Z]'))
  
  if(transposed){
    res_parse <- res_parse |>
      mutate(ind_index = col_index,
             col_index = row_index,
             row_index = ind_index,
             ind_index = NULL)
  }
  
  g1 <- res_parse |>
    mutate(row_index = fct_inorder(row_index) |> fct_rev()) |>
    ggplot(aes(col_index,row_index,fill = Cp, label = Cp)) +
    geom_tile() +
    geom_text() +
    scale_fill_gradient(low = 'orange', high = 'grey')
  
  print(g1)
  
  if(transposed){
    res_parse <- res_parse |>
      rowwise() |>
      mutate(row_index = (row_index - 1) %/% 3) |>
      ungroup()
  } else {
    res_parse <- res_parse |>
      rowwise() |>
      mutate(row_index = (utf8ToInt(row_index) + 1) %/% 3) |>
      ungroup()
  }
  
  if(is.null(group)){
    group <-
      str_c('group' ,1:length(unique(res_parse$col_index)))
  }
  
  res_parse <- res_parse |>
    distinct(col_index) |>
    mutate(group = group) |>
    right_join(res_parse)
  
  if(is.null(gene)){
    gene <-
      str_c('gene' ,1:length(unique(res_parse$row_index)))
  }
  
  res_parse |>
    distinct(row_index) |>
    arrange(row_index) |>
    mutate(gene = gene) |>
    right_join(res_parse)
}

normalize_to_refgene <- function(data, ref_gene, remove_ref = TRUE) {
  data <- data |>
    filter(gene == ref_gene) |>
    summarise(ref_cp = mean(Cp), .by = group) |>
    right_join(data) |>
    mutate(relexpr = 2^(ref_cp - Cp)) |>
    select(gene, group, relexpr)
  
  data <- data |>
    filter(gene == ref_gene) |>
    summarise(ref_exp = mean(relexpr), .by = group) |>
    right_join(data) |>
    mutate(expr = relexpr / ref_exp) |>
    select(gene, group, expr)
  
  if(remove_ref){
    data <- data |> filter(gene != ref_gene)
  }
  
  data
}

qp.data <- '00_util_scripts/qpcr/data/'

# 2405 U937 TT IgG dose ---------
dt.2405v2 <- read_csv('00_util_scripts/qpcr/data/WellResult-2024-05-07-11215.csv')

dt.2405v2 <- dt.2405v2 |>
  filter(str_detect(Cq, '^\\d')) |>
  select(Sample,Target,Cq) |>
  mutate(Cq = as.double(Cq))

dt.2405v2 |> filter(Target == 'ARG1', Sample == 'blank')
  ggplot(aes(Sample, Cq)) +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~Target)

exp.2405 <- dt.2405v2 |>
  filter(Target == 'GAPDH') |>
  summarise(ref.cq = mean(Cq), .by = Sample) |>
  left_join(dt.2405v2) |>
  mutate(rel.exp = 2^(ref.cq-Cq)) |>
  filter(!(Target == 'ARG1' & rel.exp > .01))

exp.2405 <- exp.2405 |>
  filter(Target == 'GAPDH') |>
  summarise(scale.fct = mean(rel.exp), .by = Sample) |>
  left_join(exp.2405) |>
  mutate(rel.exp = rel.exp / scale.fct)

exp.2405 |>
  filter(Target != 'GAPDH') |>
  mutate(Sample = str_replace(Sample, '-', 'ng-'),
        Sample = ifelse(str_ends(Sample, '\\d'), str_c(Sample, 'h'), Sample),
         Sample = fct_relevel(Sample, 'blank', 'IFN', 'LPS')) |>
  ggplot(aes(Sample, rel.exp)) +
  geom_jitter(height = 0, width = .1) +
  stat_summary(fun = 'mean', geom = 'crossbar', color = 'red') +
  facet_wrap(~Target, scales = 'free_y')+
  theme_pubr(x.text.angle = 45) +
  labs(title = 'TT U937 qPCR',
       y = 'Relative expression')

# 240607 U937 TT IgG time -----
parse_qpcr_cp('00_util_scripts/qpcr/data/240607-u937.txt',
              gene = c('GAPDH','IL1B','IL6','TNF','TGFB1'),
              group = c('2h','2h','4h','4h','6h','6h','8h','8h',
                        'IFN','LPS','NC')) |>
  normalize_to_refgene('GAPDH') |>
  mutate(group = fct_relevel(group, 'NC','IFN','LPS')) |>
  ggplot(aes(group, expr)) +
  geom_jitter(height = 0, width = .1) +
  stat_summary(fun = 'mean', geom = 'crossbar', color = 'red') +
  facet_wrap(~gene, scales = 'free_y')+
  theme_pubr(x.text.angle = 45) +
  labs(title = 'TT U937 qPCR',
       y = 'Relative expression')

# 240629 U937 TT IgG dose------
group.0629 <- read_tsv('00_util_scripts/qpcr/data/240629.grouping.txt',
                       col_names = 'X1')

dose29 <- parse_qpcr_cp('00_util_scripts/qpcr/data/240629-u937.txt',
              gene = c('GAPDH','IL1B','IL6','TNF','ARG1','X'),
              group = group.0629$X1)

dose.td <- dose29 |> filter(gene != 'X') |>
  normalize_to_refgene('GAPDH') |>
  mutate(dose = str_extract(group, '\\d+') |> as.integer() |> replace_na(-1),
         group = fct_reorder(group, dose) |> fct_relevel('NC','T.NC'))

dose.td |>
  filter(!str_detect(group, '^T')) |>
  ggplot(aes(group, expr)) +
  geom_jitter(height = 0, width = .3) +
  stat_summary(fun = 'mean', geom = 'crossbar', color = 'red') +
  facet_wrap(~gene, scales = 'free_y') +
  theme_pubr(x.text.angle = 90) +
  expand_limits(y = 0) +
  labs(title = 'Vector U937', x = 'IgG dose (ug/ml)', y = 'Relative foldchange')

dose.td |>
  filter(str_detect(group, '^T')) |>
  mutate(group = str_remove(group, 'T.') |>
           fct_reorder(dose) |> fct_relevel('NC')) |>
  ggplot(aes(group, expr)) +
  geom_jitter(height = 0, width = .3) +
  stat_summary(fun = 'mean', geom = 'crossbar', color = 'red') +
  facet_wrap(~gene, scales = 'free_y') +
  theme_pubr(x.text.angle = 90) +
  expand_limits(y = 0) +
  labs(title = 'FCGR2B-232 TT U937',
       x = 'IgG dose (ug/ml)', y = 'Relative foldchange')

dose.td |>
  mutate(genotype = ifelse(str_detect(group, 'T'), 'TT', 'Vector')) |>
  ggplot(aes(group, expr)) +
  geom_jitter(height = 0, width = .3) +
  stat_summary(aes(color = genotype),fun = 'mean', geom = 'crossbar') +
  facet_wrap(~gene, scales = 'free_y') +
  theme_pubr(x.text.angle = 90) +
  expand_limits(y = 0) +
  labs(title = 'Vector & FCGR2B-232TT U937',
       x = 'IgG dose (ug/ml)', y = 'Relative foldchange')

# 250429 qPCR ------
group0429 <- read_csv('00_util_scripts/qpcr/data/250429.group.txt') |>
  pull(group)

cp0429 <-
parse_qpcr_cp('00_util_scripts/qpcr/data/20250429a431.txt',
              gene = c('GAPDH','CXCL8','IL6','CCL20'), group = group0429)

cp0429 |>
  normalize_to_refgene(ref_gene = 'GAPDH') |>
  ggplot(aes(expr, group)) +
  stat_summary(geom = 'col', fun = mean) +
  ggbeeswarm::geom_beeswarm() +
  facet_wrap(~gene)

cp0429 |>
  filter(group == 'DMSO')

# qPCR visualization ------
data <- read_tsv('data/other/220625-qP-RelaExp.txt')

data %>%
  filter(Time == '2') %>%
  ggplot(aes(x = Substance, y = RelativeExp)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  labs(title = 'IL-1B 2h')

data %>%
  filter(Time == '6') %>%
  ggplot(aes(x = Substance, y = RelativeExp)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  coord_flip()+
  labs(title = 'IL-1B 6h')

data %>%
  filter(Gene == 'IL-10' & Group == '6h') %>%
  ggplot(aes(x = Name, y = fc)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  labs(title = 'IL-10 6h')

data %>%
  filter(Gene == 'IL-1B' & Group == '2h') %>%
  ggplot(aes(x = Name, y = fc)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  labs(title = 'IL-1B 2h')

data %>%
  filter(Gene == 'IL-1B' & Group == '6h') %>%
  ggplot(aes(x = Name, y = fc)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  labs(title = 'IL-1B 6h')

data %>%
  filter(Gene == 'IL-1B') %>%
  ggplot() +
  geom_dotplot(binaxis = 'y', stackdir = 'center', aes(x = Name, y = fc, fill = Group)) +
  labs(title = 'IL-1B')

data %>%
  filter(Gene == 'TNFA' & Group == '2h') %>%
  ggplot(aes(x = Name, y = fc)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  labs(title = 'TNFA 2h')

data %>%
  filter(Gene == 'TNFA' & Group == '6h') %>%
  ggplot(aes(x = Name, y = fc)) +
  geom_dotplot(binaxis = 'y', stackdir = 'center') +
  stat_summary(geom = "crossbar", fun = mean, colour = "red") +
  ylim(0,8)+
  labs(title = 'TNFA 6h')

data %>%
  filter(Gene == 'TNFA' & Group == '6h') %>%
  ggplot() +
  geom_boxplot(aes(x = Name, y = fc)) +
  labs(title = 'TNFA 6h')

## 230727 BMDM qPCR ------------
bmdm727 <- read_delim('00_util_scripts/data/0727-qpcr.txt', skip = 1)

td_bmdm727 <- bmdm727 |>
  select(Pos, Cp) |>
  mutate(col_index = str_extract(Pos, '\\d+') |> as.numeric(),
         row_index = str_extract(Pos, '[A-Z]'),
         Cp = replace_na(Cp, 36),
         gene = case_match(row_index,
                           c('A','B','C') ~ 'Gapdh',
                           c('D','E','F') ~ 'Il-1b',
                           c('G','H','I') ~ 'Il-6',
                           c('J','K','L') ~ 'Tnf',
                           c('M','N','O') ~ 'Tgfb1'),
         sample = case_match(col_index,
                             c(1,3,5) ~ 'IgG',
                             c(7,9,11) ~ 'PBS',
                             c(13,15,17) ~ 'Soluble IgG',
                             c(19,21,23) ~ 'Medium'
         ),
         genotype = case_match(col_index,
                               c(1,7,13,19) ~ 'II',
                               c(3,9,15,21) ~ 'IT',
                               c(5,11,17,23) ~ 'TT')) |>
  filter(!is.na(sample) & row_index != 'P')

td_bmdm727 |>
  ggplot(aes(as.numeric(col_index),row_index,fill = Cp, label = Cp)) +
  geom_raster() +
  geom_text() +
  scale_fill_gradient(low = 'orange', high = 'grey')


## 230901 BMDM qPCR -----
bmdm <- read_delim('00_util_scripts/data/20230901-bmdm-qpcr.txt', skip = 1)

td_bmdm <- bmdm |>
  select(Pos, Cp) |>
  mutate(col_index = str_extract(Pos, '\\d+'),
         row_index = str_extract(Pos, '[A-Z]'),
         Cp = replace_na(Cp, 36),
         gene = case_match(row_index,
                           c('A','B','C') ~ 'Gapdh',
                           c('D','E','F') ~ 'Il-1b',
                           c('G','H','I') ~ 'Il-6',
                           c('J','K','L') ~ 'Tnf',
                           c('M','N','O') ~ 'Tgfb1'),
         sample = case_match(col_index,
                             c('1','3','5') ~ 'IgG',
                             c('7','9','11') ~ 'PBS',
                             c('13','15','17') ~ 'Soluble IgG',
         ),
         genotype = case_match(col_index,
                               c('1','7','13','19') ~ 'II',
                               c('3','9','15','21') ~ 'IT',
                               c('5','11','17','23') ~ 'TT')) |>
  filter(!is.na(sample) & row_index != 'P')

td_bmdm |>
  ggplot(aes(as.numeric(col_index),row_index,fill = Cp, label = Cp)) +
  geom_raster() +
  geom_text() +
  scale_fill_gradient(low = 'orange', high = 'grey')

## PBS group Gapdh too low, not usable in this experimemt
td_bmdm <- td_bmdm |>
  filter(sample != 'PBS')

res_bmdm <- td_bmdm |>
  filter(gene == 'Gapdh') |>
  summarise(ref_cp = mean(Cp), .by = col_index) |>
  right_join(td_bmdm) |>
  mutate(foldchange = 2^(ref_cp - Cp)) |>
  select(gene, sample, genotype, foldchange) |>
  filter(gene != 'Gapdh')

res_bmdm

res_bmdm |>
  filter(sample == 'IgG') |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr()

res_bmdm |>
  filter(sample == 'Soluble IgG') |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr()

ctl_bmdm <- res_bmdm |>
  summarise(mean_fc = mean(foldchange), .by = c(gene, sample, genotype)) |>
  filter(sample == 'PBS') |>
  rename(ctrl.fc = mean_fc) |>
  select(-sample) |>
  right_join(res_bmdm)

ctl_bmdm |>
  filter(sample == 'IgG' & gene != 'Tgfb1') |>
  mutate(foldchange = foldchange/ctrl.fc) |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  geom_hline(yintercept = 2, linetype = 'dashed') +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr() +
  labs(title = 'BMDM cytokine expression after 24h IgG stimulation',
       y = 'Fold change') +
  stat_compare_means(comparisons = list(c('TT','IT'),
                                        c('TT','II')),
                     method = 't.test') +
  scale_color_manual(values = c('green3','blue','red'))

g1 <- ctl_bmdm |>
  filter(sample == 'IgG' & gene == 'Il-1b') |>
  mutate(foldchange = foldchange/ctrl.fc) |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  geom_hline(yintercept = 2, linetype = 'dashed') +
  theme_pubr() +
  labs(y = 'IL-1B Fold change') +
  stat_compare_means(comparisons = list(c('TT','IT'),
                                        c('TT','II')),
                     method = 't.test') +
  scale_color_manual(values = c('green3','blue','red'))

g2 <- ctl_bmdm |>
  filter(sample == 'IgG' & gene == 'Il-6') |>
  mutate(foldchange = foldchange/ctrl.fc) |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  geom_hline(yintercept = 2, linetype = 'dashed') +
  theme_pubr() +
  labs(y = 'IL-6 Fold change') +
  stat_compare_means(comparisons = list(c('TT','IT'),
                                        c('TT','II')),
                     method = 't.test') +
  scale_color_manual(values = c('green3','blue','red'))

g3 <- ctl_bmdm |>
  filter(sample == 'IgG' & gene == 'Tnf') |>
  mutate(foldchange = foldchange/ctrl.fc) |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  geom_hline(yintercept = 2, linetype = 'dashed') +
  theme_pubr() +
  labs(y = 'TNFa Fold change') +
  stat_compare_means(comparisons = list(c('TT','IT'),
                                        c('TT','II')),
                     method = 't.test') +
  scale_color_manual(values = c('green3','blue','red'))

g1 + g2 + g3 + patchwork::plot_layout(guides = 'collect') +
  patchwork::plot_annotation(theme = theme(legend.position = 'top'),
                             title = 'BMDM cytokine expression after 24h IgG stimulation',)

## 230902 BMDM qPCR ------
bmdm92 <- read_delim('00_util_scripts/data/20230902-bmdm-qpcr.txt', skip = 1)

td_bmdm92 <- bmdm92 |>
  select(Pos, Cp) |>
  mutate(col_index = str_extract(Pos, '\\d+') |> as.numeric(),
         row_index = str_extract(Pos, '[A-Z]'),
         gene = case_match(row_index,
                           c('A','B','C') ~ 'Gapdh',
                           c('D','E','F') ~ 'Il-1b',
                           c('G','H','I') ~ 'Il-6',
                           c('J','K','L') ~ 'Tnf',
                           c('M','N','O') ~ 'Tgfb1'),
         sample = case_match(col_index,
                             c(1,3,5) ~ 'IgG',
                             c(7,9,11) ~ 'PBS',
         ),
         genotype = case_match(col_index,
                               c(1,7) ~ 'II',
                               c(3,9) ~ 'IT',
                               c(5,11) ~ 'TT')) |>
  filter(!is.na(gene) & !is.na(sample))

td_bmdm92 |>
  ggplot(aes(col_index,row_index,fill = Cp, label = Cp)) +
  geom_raster() +
  geom_text() +
  scale_fill_gradient(low = 'orange', high = 'grey')

## although some gene is low, but all Gapdh is high
td_bmdm92 <- td_bmdm92 |>
  mutate(Cp = replace_na(Cp, 36))

res_bmdm92 <- td_bmdm92 |>
  filter(gene == 'Gapdh') |>
  summarise(ref_cp = mean(Cp), .by = col_index) |>
  right_join(td_bmdm92) |>
  mutate(foldchange = 2^(ref_cp - Cp)) |>
  select(gene, sample, genotype, foldchange) |>
  filter(gene != 'Gapdh')

res_bmdm92

res_bmdm92 |>
  filter(sample == 'IgG') |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr()

res_bmdm92 |>
  filter(sample == 'PBS') |>
  ggplot(aes(genotype, foldchange, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr()

relexpr_92 <- res_bmdm92 |>
  filter(sample == 'PBS') |>
  group_by(gene, genotype) |>
  summarise(ref_fc = mean(foldchange)) |>
  right_join(res_bmdm92) |>
  mutate(expr = foldchange / ref_fc) |>
  select(gene, sample, genotype, expr) |>
  filter(sample != 'PBS')

relexpr_92 |>
  filter(gene == 'Tnf') |>
  ggplot(aes(genotype, expr, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  labs(y = 'TNFa foldchange relative to PBS', title = 'TNFa') +
  theme_pubr() +
  stat_compare_means(comparisons = list(c('TT','IT'),c('II','TT')),
                     method = 't.test')

g1 <- last_plot()

hybrid91.92 <- res_bmdm92 |>
  filter(sample == 'PBS') |>
  group_by(gene, genotype) |>
  summarise(ref_fc = mean(foldchange)) |>
  right_join(res_bmdm) |>
  mutate(expr = foldchange / ref_fc) |>
  select(gene, sample, genotype, expr) |>
  filter(sample != 'PBS')

hybrid91.92 |>
  filter(sample == 'Soluble IgG' & gene == 'Il-1b') |>
  ggplot(aes(genotype, expr, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  labs(y = 'Il-1b foldchange relative to PBS', title = 'IL-1b') +
  theme_pubr() +
  stat_compare_means(comparisons = list(c('TT','IT'),c('II','TT')),
                     method = 't.test')

g2 <- last_plot()

hybrid91.92 |>
  filter(sample == 'IgG' & gene == 'Tgfb1') |>
  ggplot(aes(genotype, expr, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  labs(y = 'TGFb1 foldchange relative to PBS', title = 'TGFb1') +
  theme_pubr() +
  stat_compare_means(comparisons = list(c('TT','IT'),c('II','TT')),
                     method = 't.test')

g3 <- last_plot()

g1 + g2 + g3 + plot_layout(guides = 'collect') +
  plot_annotation(theme = theme(legend.position = 'top'),
                  title = 'BMDM cytokine expression after 24h IgG stimulation',)

## 230920 U937 FCGR2 expr ------
u937 <- parse_qpcr_cp('00_util_scripts/data/20230920-u937-fcgr2.txt',
                      gene = c('GAPDH','FCGR2A','FCGR2B','FCGR2C'),
                      group = c('II','vector','TT')) |> 
  normalize_to_refgene('GAPDH')

u937 |>
  ggplot(aes(group, relexpr, color = group)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene) +
  theme_pubr()

## 230921 BMDM IC stimulation ------
genotype_group <- expand_grid(a = c('II','IT','TT'), b = c('IC1','IC2','IgG','PBS')) |>
  unite(c, a:b)

bmdmic <- parse_qpcr_cp('00_util_scripts/data/20230921-bmdm-ic.txt',
                        gene = c('GAPDH','IL1B','IL6','TNFA','TGFB'),
                        group = genotype_group$c) |> 
  normalize_to_refgene('GAPDH')

bmdmic <- bmdmic |>
  separate(group,into = c('genotype','treatment'))

bmdmic_res <- bmdmic |>
  filter(treatment == 'PBS') |>
  group_by(gene, genotype) |>
  summarise(ref_fc = mean(relexpr)) |>
  right_join(bmdmic) |>
  mutate(expr = relexpr / ref_fc) |>
  select(gene, treatment, genotype, expr)

bmdmic_res |>
  filter(treatment == 'IgG' & gene != 'GAPDH') |>
  ggplot(aes(genotype, expr, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr() +
  stat_compare_means(comparisons = list(c('TT','IT'),c('II','TT')),
                     method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0,.1))) +
  labs(y = 'Expression fold change to PBS', title = 'IgG stimulation of BMDM')

bmdmic_res |>
  filter(str_detect(treatment, 'IC1') & gene != 'GAPDH') |>
  ggplot(aes(genotype, expr, color = genotype)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr() +
  stat_compare_means(comparisons = list(c('TT','IT'),c('II','TT')),
                     method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0,.1))) +
  labs(y = 'Expression fold change to PBS', title = 'IgG-beads stimulation of BMDM')

## qPCR of ycj U937 ----
ycjcd <- read_delim('00_util_scripts/data/u937-ycj-qp.txt',
                    col_names = c('group','rel'), col_types = 'cd')

ycjcd |>
  ggplot(aes(group, rel, color = group)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'Expression of FCGR2B mRNA relative to vector',
       y = 'Relative expression')

## 240407 -----------
group0407 <- c('blank','IFNg+crosslink','blank','IFNg+crosslink','LPS','IFNg+LPS','LPS','IFNg+LPS',
               'IFNg','IFNg','IFNg+IgG','IFNg+IgG')

u937tt <- parse_qpcr_cp('00_util_scripts/data/240407-qpcr-u937.txt', transposed = T,
                        gene = c('GAPDH','IL1B','IL6','ARG1','NOS2','MRC1','VEGFA'),
                        group = group0407)

u937tt |>
  normalize_to_refgene('GAPDH',remove_ref = T) |>
  filter(gene != 'NOS2') |>
  ggplot(aes(group, expr, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr(x.text.angle = 45, legend = 'none') +
  scale_color_manual(values = c('black','blue','orange','red','purple','green3')) +
  labs(y = 'Expression relative to GAPDH') +
  facet_wrap(~gene, scales = 'free_y')

plotly::ggplotly()

## 231105 --------
u937cyt <- read_delim('00_util_scripts/data/231008-u937-igg.txt')

u937prs <- u937cyt |>
  fill(gene) |>
  pivot_longer(2:last_col(), names_to = 'group') |>
  mutate(ctl = ifelse(str_detect(group, 'without'), 'ctl', 'exp'),
         group = str_extract(group, 'U937|II|TT'))

u937rel <- u937prs |>
  filter(gene == 'GAPDH') |>
  group_by(group, ctl) |>
  summarise(ref_gene = mean(value)) |>
  right_join(u937prs) |>
  mutate(relexpr = 2^(ref_gene-value))

u937rel |>
  filter(ctl == 'ctl') |>
  group_by(group,gene) |>
  summarise(ctl_expr = mean(relexpr)) |>
  right_join(u937rel) |>
  filter(ctl == 'exp' & gene != 'GAPDH') |>
  ungroup() |>
  mutate(foldchange = relexpr / ctl_expr,
         group = fct_relevel(group, 'U937', 'II')) |>
  ggplot(aes(group, foldchange, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr() +
  labs_pubr() +
  labs(x = 'group', y = 'relative expression',
       title = '10ug/ml IgG stimulation') +
  stat_compare_means(comparisons = list(c('II','TT'),
                                        c('U937','TT')),
                     method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0, .2))) +
  scale_color_manual(values = c('black','blue','red'))

u937rel |>
  filter(ctl == 'ctl') |>
  group_by(group,gene) |>
  summarise(ctl_expr = mean(relexpr)) |>
  right_join(u937rel) |>
  filter(ctl == 'exp' & gene != 'GAPDH' & group != 'U937') |>
  ungroup() |>
  mutate(foldchange = relexpr / ctl_expr,
         group = fct_relevel(group, 'II')) |>
  ggplot(aes(group, foldchange, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene, scales = 'free') +
  theme_pubr() +
  labs_pubr() +
  labs(x = 'group', y = 'relative expression',
       title = '10ug/ml IC stimulation') +
  stat_compare_means(comparisons = list(c('II','TT')),
                     method = 't.test') +
  scale_y_continuous(expand = expansion(mult = c(0, .2))) +
  scale_color_manual(values = c('blue','red'))